#include "Material.h"

Material::Material()
{
    _groups = 0;
    _sigma_t = NULL;
    _sigma_a = NULL;
    _sigma_s = NULL;
    _sigma_f = NULL;
    _nu = NULL;
    _nu_sigma_f = NULL;
    _chi = NULL;
    // _fmtin = MCNP;
    // _fmtout = MCNP;

    _id = 0;
    _up = 0;
    _down = 0;
    _density = 0.0;
}


Material::~Material()
{
    if (_sigma_t != NULL) delete [] _sigma_t;
    if (_sigma_a != NULL) delete [] _sigma_a;
    if (_sigma_s != NULL) delete [] _sigma_s;
    if (_sigma_f != NULL) delete [] _sigma_f;
    if (_nu_sigma_f != NULL) delete [] _nu_sigma_f;
    if (_nu != NULL) delete [] _nu;
    if (_chi != NULL) delete [] _chi;
    if (_energy_center != NULL) delete [] _energy_center;
    if (_energy_width != NULL) delete [] _energy_width;
}


// void Material::setInFMT(datfmt in) { _fmtin = in; }
// void Material::setOutFMT(datfmt out) { _fmtout = out; }
void Material::setName(std::string name)
{
    int length = name.length();
    if (length <= 0)
    {
        std::cerr << "error: setName(), string is empty." << std::endl;
        return;
    }
    _name = new char [length+1];
    strcpy(_name, name.c_str());
    _name[length] = '\0';
}


void Material::setGroups(std::vector<int> v)
{
    if (v.size() <= 0)
    {
        std::cerr << "error: setGroups(), vector is empty." << std::endl;
        return;
    }

    _groups = v[0];
}


void Material::setID(std::vector<int> v)
{
    if (v.size() <= 0)
    {
        std::cerr << "error: setID(), vector is empty." << std::endl;
        return;
    }

    _id = v[0];
}


void Material::setUp(std::vector<int> v)
{
    if (v.size() <= 0)
    {
        std::cerr << "error: setUp(), vector is empty." << std::endl;
        return;
    }

    _up = v[0];
}


void Material::setDown(std::vector<int> v)
{
    if (v.size() <= 0)
    {
        std::cerr << "error: setDown(), vector is empty." << std::endl;
        return;
    }

    _down = v[0];
}


void Material::setDensity(std::vector<double> v)
{
    if (v.size() <= 0)
    {
        std::cerr << "error: setDensity(), vector is empty." << std::endl;
        return;
    }

    _density = v[0];
}


void Material::setSigma_t(std::vector<double> v)
{
    if (v.size() <= 0)
    {
        std::cerr << "error: setSimga_t(), vector is empty." << std::endl;
        return;
    }

    if (v.size() == _groups)
    {
        for (int i=0; i<_groups; i++)
        {
            _sigma_t[i] = v[i];
        }
    }
    else
    {
        std::cerr << "error: setSigma_t(), vector size does not match energy groups." << std::endl;
        return;
    }
}


void Material::setSigma_f(std::vector<double> v)
{
    if (v.size() <= 0)
    {
        std::cerr << "error: setSimga_f(), vector is empty." << std::endl;
        return;
    }

    if (v.size() == _groups)
    {
        for (int i=0; i<_groups; i++)
        {
            _sigma_f[i] = v[i];
        }
    }
    else
    {
        std::cerr << "error: setSigma_f(), vector size does not match energy groups." << std::endl;
        return;
    }
}


void Material::setSigma_a(std::vector<double> v)
{
    if (v.size() <= 0)
    {
        std::cerr << "error: setSimga_a(), vector is empty." << std::endl;
        return;
    }

    if (v.size() == _groups)
    {
        for (int i=0; i<_groups; i++)
        {
            _sigma_a[i] = v[i];
        }
    }
    else
    {
        std::cerr << "error: setSigma_a(), vector size does not match energy groups." << std::endl;
        return;
    }
}


void Material::setSigma_s(std::vector<double> v, std::string string)
{
    // std::cout << v.size() << std::endl;
    // std::cout << _up << std::endl;
    // std::cout << _down << std::endl;
    // std::cout << _groups*(1+_up+_down)-_down*(_down+1)/2-_up*(_up+1)/2 << std::endl;
    if (v.size() <= 0)
    {
        std::cerr << "error: setSimga_s(), vector is empty." << std::endl;
        return;
    }

    if (v.size() == _groups*_groups)
    {
        for (int i=0; i<_groups; i++)
        {
            for(int j=0; j<_groups; j++)
            {
                if (string == "upper")
                {
                    _sigma_s[j*_groups+i] = v[i*_groups+j]; // upper
                }
                else if (string == "lower")
                {
                    _sigma_s[j*_groups+i] = v[j*_groups+i]; // lower
                }
                else
                {
                    std::cout << "scattering_matrix error: " << string << std::endl;
                }
            }
        }
    }
    else if (v.size() == _groups*(1+_up+_down)-_down*(_down+1)/2-_up*(_up+1)/2)
    {
        int counter = 0;
        for (int i=0; i<_groups; i++)
            for (int j=0; j<_groups; j++)
            {
                if ((j+1) >= max(1, i+1-_up) && (j+1) <= min(_groups, i+1+_down))
                {
                    _sigma_s[j*_groups+i] = v[counter++];
                }
            }
    }
    else
    {
        std::cerr << "error: setSigma_s(), vector size does not match energy groups." << std::endl;
        return;
    }
}


void Material::setChi(std::vector<double> v)
{
    if (v.size() <= 0)
    {
        std::cerr << "error: setChi(), vector is empty." << std::endl;
        return;
    }

    if (v.size() == _groups)
    {
        for (int i=0; i<_groups; i++)
        {
            _chi[i] = v[i];
        }
    }
    else
    {
        std::cerr << "error: setChi(), vector size does not match energy groups." << std::endl;
        return;
    }
}


void Material::setNu(std::vector<double> v)
{
    if (v.size() <= 0)
    {
        std::cerr << "error: setNu(), vector is empty." << std::endl;
        return;
    }

    if (v.size() == _groups)
    {
        for (int i=0; i<_groups; i++)
        {
            _nu[i] = v[i];
        }
    }
    else
    {
        std::cerr << "error: setNu(), vector size does not match energy groups." << std::endl;
        return;
    }
}


void Material::setEnergyCenter(std::vector<double> v)
{
    if (v.size() <= 0)
    {
        std::cerr << "error: setEnergyCenter(), vector is empty." << std::endl;
        return;
    }

    if (v.size() == _groups)
    {
        for (int i=0; i<_groups; i++)
        {
            _energy_center[i] = v[i];
        }
    }
    else
    {
        std::cerr << "error: setEnergyCenter(), vector size does not match energy groups." << std::endl;
        return;
    }
}


void Material::setEnergyWidth(std::vector<double> v)
{
    if (v.size() <= 0)
    {
        std::cerr << "error: setEnergyWidth(), vector is empty." << std::endl;
        return;
    }

    if (v.size() == _groups)
    {
        for (int i=0; i<_groups; i++)
        {
            _energy_width[i] = v[i];
        }
    }
    else
    {
        std::cerr << "error: setEnergyWidth(), vector size does not match energy groups." << std::endl;
        return;
    }
}


void Material::initialize()
{
    if (_groups <= 0)
        std::cout << "Error: Can NOT allocate memory for the Material, "
            << "cause the number of groups is no more than 0.\n";

    _sigma_t = new double [_groups];
    _sigma_a = new double [_groups];
    _sigma_s = new double [_groups*_groups];
    _sigma_f = new double [_groups];
    _nu_sigma_f = new double [_groups];
    _nu = new double [_groups];
    _chi = new double [_groups];
    _energy_center = new double [_groups];
    _energy_width = new double [_groups];
}


void Material::convert2Macro()
{
    /* compute macro crossection */
    if (_density == 0.0)
    {
        for (int i=0; i<_groups; i++)
        {
            _nu_sigma_f[i] = _nu[i]*_sigma_f[i];
        }
    }
    else
    {
        for (int i=0; i<_groups; i++)
        {
            _sigma_t[i] = _density*_sigma_t[i];
            _sigma_f[i] = _density*_sigma_f[i];
            _sigma_a[i] = _density*_sigma_a[i];
            _nu_sigma_f[i] = _density*_nu[i]*_sigma_f[i];
        }
        for (int i=0; i<_groups*_groups; i++)
        {
            _sigma_s[i] = _density*_sigma_s[i];
        }
    }
}


void Material::print()
{
    std::cout << "Material name: " << _name << ", Material id: " << _id;
    std::cout << "\nGroup number: " << _groups;
    std::cout << "\nCross Sections:\n";
    std::cout << "Sigma_t: ";
    for (int i=0; i<_groups; i++)
        std::cout << _sigma_t[i] << ", ";
    std::cout << "\nSigma_a: ";
    for (int i=0; i<_groups; i++)
        std::cout << _sigma_a[i] << ", ";
    std::cout << "\nSigma_s:\n";
    for (int i=0; i<_groups; i++)
    {
        std::cout << "         ";
        for (int j=0; j<_groups; j++)
            std::cout << _sigma_s[j*_groups+i] << ", ";
        std::cout << "\n";
    }
    std::cout << "Sigma_f: ";
    for (int i=0; i<_groups; i++)
        std::cout << _sigma_f[i] << ", ";

    std::cout << "\nNu: ";
    for (int i=0; i<_groups; i++)
        std::cout << _nu[i] << ", ";

    std::cout << "\nChi: ";
    for (int i=0; i<_groups; i++)
        std::cout << _chi[i] << ", ";
    std::cout << std::endl;
}


void Material::write(std::ofstream & fout, datfmt fmt)
{
    switch(fmt)
    {
        /* write data to XML file */
        case XML:
        {
            /* write the material id and name */
            fout << "<material id=\"" << _id << "\" name=\"" << _name << "\">\n";

            /* wirte the cross section and number of groups */
            fout << "<cross_section group=\"" << _groups << "\">\n";

            /* write sigma_t */
            fout << "<sigma_t>";
            for (int i=0; i<_groups; i++)
            {
                if (i == 0)
                    fout << std::setiosflags(std::ios::scientific) << _sigma_t[i];
                else
                    fout << std::setw(20) << std::setiosflags(std::ios::scientific)
                    << _sigma_t[i];
            }
            fout << "</sigma_t>\n";

            /* write sigma_s */
            fout << "<sigma_s>";
            for (int i=0; i<_groups; i++)
                for (int j=0; j<_groups; j++)
                {
                    if (j*_groups+i == 0)
                        fout << std::setiosflags(std::ios::scientific) << _sigma_s[j*_groups+i];
                    else
                        fout << std::setw(20) << std::setiosflags(std::ios::scientific)
                        << _sigma_s[j*_groups+i];
                }
            fout << "</sigma_s>\n";

            /* write sigma_f */
            fout << "<sigma_f>";
            for (int i=0; i<_groups; i++)
            {
                if (i == 0)
                    fout << std::setiosflags(std::ios::scientific) << _sigma_f[i];
                else
                    fout << std::setw(20) << std::setiosflags(std::ios::scientific)
                    << _sigma_f[i];
            }
            fout << "</sigma_f>\n";

            /* write nu_sigma_f */
            fout << "<nu_sigma_f>";
            for (int i=0; i<_groups; i++)
            {
                if (i == 0)
                    fout << std::setiosflags(std::ios::scientific) << _nu_sigma_f[i];
                else
                    fout << std::setw(20) << std::setiosflags(std::ios::scientific)
                    << _nu_sigma_f[i];
            }
            fout << "</nu_sigma_f>\n";

            /* write chi */
            fout << "<chi>";
            for (int i=0; i<_groups; i++)
            {
                if (i == 0)
                    fout << std::setiosflags(std::ios::scientific) << _chi[i];
                else
                    fout << std::setw(20) << std::setiosflags(std::ios::scientific)
                    << _chi[i];
            }
            fout << "</chi>\n";
            fout << "</cross_section>\n"; /* end of cross section */
            fout << "</material>\n"; /* end of material */
            break;
        }
        /* wirte data to MCNP file */
        case MCNP:
        {
            /* nxs, jxs array, set the array according MCNP manual */
            int nxs[17], jxs[33];
            for (int i=0; i<33; i++)
            {
                if (i < 17)
                    nxs[i] = 0;
                jxs[i] = 0;
            }

            /* number of groups */
            nxs[5] = _groups;

            /* upscater groups */
            nxs[6] = _up;

            /* downscatter groups */
            nxs[7] = _down;

            /* number of nubars given */
            nxs[10] = 1;

            /* identifier for incident particle
            IPT = 1 for neutrons
            IPT = 2 for photons
            IPT = 0 for other particles(temporary)
             */
            nxs[12] = 1;

            /* length of second block of data */
            nxs[1] = 2*nxs[5] + nxs[5] + nxs[5]*nxs[10] + nxs[5] + nxs[5] + nxs[5] +
                nxs[5]*(1+nxs[7]+nxs[6])-nxs[7]*(nxs[7]+1)/2-nxs[6]*(nxs[6]+1)/2+ 3;

            /* ZA = 1000 * Z + A */
            nxs[2] = _id;

            /* location of incident particle group structure = 1 */
            jxs[1] = 1;

            /* lcation of total xs */
            jxs[2] = jxs[1] + 2*nxs[5];

            /* location of fission xs */
            jxs[3] = jxs[2] + nxs[5];

            /* location of nubar data */
            jxs[4] = jxs[3] + nxs[5]*nxs[10];

            /* location of fission chi data */
            jxs[5] = jxs[4] + nxs[5];

            /* location of absorption xs */
            jxs[6] = jxs[5] + nxs[5];

            /* location of P0 locators */
            jxs[13] = jxs[5] + nxs[5] + nxs[5]
                      + nxs[5]*(1+nxs[7]+nxs[6])-nxs[7]*(nxs[7]+1)/2-nxs[6]*(nxs[6]+1)/2;

            /* location of XPN locators */
            jxs[16] = jxs[13] + nxs[8] + 1;

            /* location of PN locators */
            jxs[17] = jxs[16] + nxs[8] + 1;

            /* get time */
            time_t rawtime;
            struct tm* timeinfo;
            char   buffer[32];
            time(&rawtime);
            timeinfo = localtime(&rawtime);
            strftime(buffer, 32, "%D", timeinfo);

            /* write data to MCNP data file */
            fout << std::setw(10) << _name;
            fout << std::setw(12) << std::fixed << std::setprecision(6) << _density;
            fout << std::setw(12) << std::scientific << std::setprecision(6) << 0.0
                 << std::setw(10) << buffer << "\n";


            fout << "Material "<< std::setw(5) << _id << "\n";

            for (int i=1; i<17; i++)
            {
                fout << std::setw(7) << 0
                     << std::setw(11) << std::fixed << std::setprecision(6) << 0.0;
                if (i % 4 == 0)
                    fout << "\n";
            }
            fout << std::resetiosflags(std::ios::fixed);

            /* write the NXS arrays */
            for (int i=1; i<17; i++)
            {
                fout << std::setw(9) << nxs[i];
                if (i % 8 == 0)
                    fout << "\n";
            }

            /* write the JXS arrays */
            for (int i=1; i<33; i++)
            {
                fout << std::setw(9) << jxs[i];
                if (i % 8 == 0)
                    fout << "\n";
            }

            /* write the second data block */
            double* block= new double [nxs[1]-3];

            /* enery center and width */
            for (int i=jxs[1]; i<jxs[2]; i++)
            {
                if (i-jxs[1] < _groups)
                    block[i] = _energy_center[i-jxs[1]];
                else
                    block[i] = _energy_width[i-jxs[1]-nxs[5]];
            }

            /* sigma_t */
            for (int i=jxs[2]; i<jxs[3]; i++)
            {
                block[i] = _sigma_t[i-jxs[2]];
            }

            /* sigma_f */
            for (int i=jxs[3]; i<jxs[4]; i++)
            {
                block[i] = _sigma_f[i-jxs[3]];
            }

            /* nu */
            for (int i=jxs[4]; i<jxs[5]; i++)
            {
                block[i] = _nu[i-jxs[4]];
            }

            /* chi */
            for (int i=jxs[5]; i<jxs[6]; i++)
            {
                block[i] = _chi[i-jxs[5]];
            }

            /* sigma_a */
            for (int i=jxs[6]; i<jxs[6]+nxs[5]; i++)
            {
                block[i] = _sigma_a[i-jxs[6]];
            }

            /* sigma_s */
            int counter = jxs[6] + nxs[5];
            for (int i=0; i<nxs[5]; i++)
            {
                for (int j=0; j<nxs[5]; j++)
                {
                    if ((j+1) >= max(1, i+1-nxs[6]) && (j+1) <= min(nxs[5], i+1+nxs[7]))
                    {
                        block[counter++] = _sigma_s[j*nxs[5]+i];
                    }
                }
            }

            for (int i=1; i<=nxs[1]-3; i++)
            {
                fout << std::setw(20) << std::setprecision(12) << std::scientific << block[i];

                if (i % 4 == 0)
                    fout << "\n";
            }

            for (int i=0; i<3; i++)
            {
                if (i == 0)
                    fout << std::setw(20) << jxs[6]+nxs[5];
                else
                    fout << std::setw(20) << 0;
                if ((i+1+nxs[1]-3) % 4 == 0)
                    fout << "\n";
            }

            if (nxs[1] % 4 != 0)
                fout << std::endl;

            // for (int i=jxs[6]+nxs[5]; i<jxs[13]; i++)
            // {
            //     block[i] = _sigma_s[((i-jxs[6])%_groups)*_groups+((i-jxs[6])/_groups)];
            // }
            delete [] block;
            break;
        }
    }
}


void Material::write(char* filename)
{
    if (filename != NULL && strlen(filename) > 0)
    {
        FILE* fmocp = fopen(filename, "a");

        // fprintf(fmocp, "[Materials]\n");

        fprintf(fmocp, "\t[./%dm]\n", _id);

        fprintf(fmocp, "\t\tid = %d\n", _id);
        fprintf(fmocp, "\t\tnumber_groups = %d\n", _groups);

        fprintf(fmocp, "\t\tsigma_t =");
        for(int i=0; i<_groups; i++)
            fprintf(fmocp, " %f", _sigma_t[i]);

        fprintf(fmocp, "\n\t\tnu_sigma_f =");
        for(int i=0; i<_groups; i++)
            fprintf(fmocp, " %f", _sigma_f[i] * _nu[i]);

        fprintf(fmocp, "\n\t\tchi =");
        for(int i=0; i<_groups; i++)
            fprintf(fmocp, " %f", _chi[i]);

        fprintf(fmocp, "\n\t\tsigma_s =");
        for (int i=0; i<_groups; i++)
        {
            for(int j=0; j<_groups; j++)
            {
                fprintf(fmocp, " %f", _sigma_s[j*_groups+i]);
            }
        }

        fprintf(fmocp, "\n\t[../]\n");

        // fprintf(fmocp, "[]\n", );
        fclose(fmocp);
    }
    else
    {
        printf("error: please check the MOCP type material file.\n");
    }
}




int max(int a, int b)
{
    if (a < b)
        return b;
    else
        return a;
}


int min(int a, int b)
{
    if (a < b)
        return a;
    else
        return b;
}
